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Abstract 

The evolution of large amplitude Tollmien-Schlichting waves in a supersonic 
boundary layer is investigated. Disturbances which have their wavenumber and 
frequency slowly varying in time and space are described using a phase equation 
type of approach. Unlike the incompressible case we find that the initial bifurca- 
tion to a finite amplitude Tollmien-Schlichting wave is subcritical for most Mach 
numbers. In fact the bifurcation is only supercritical for a small range of Mach 
numbers and even then for only a finite range of wave propagation angles. The 
modulational instability of large amplitude wavetrains is considered and is shown 
o je governed by an equation similar to Burgers equation but with the viscous 
erm replaced by a fractional derivative. A numerical investigation of the solution 
of this equation is described. It is shown that uniform wavetrains are unstable. 
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1 Introduction and triple-deck theory for compress- 
ible flows 

Consider the flow of a viscous compressible fluid which has pressure, density and speed 
and sufficiently far away from the semi-infinite plate defined by y = 0 ,r > 0 
wit i respej t to a Cartesian coordinate system made dimensionless using the reference 
engtl, i. If the viscosity in the free stream is < then we define a Reynolds number R 

Oy WI ltlllg 

[{ — VqqL 

V 'oo ' (b 

and throughout our investigation we assume R » 1. We scale the velocity and tem- 
perature of the flow near the wall on their free stream values and we assume a linear 
viscosity temperature law which we write in the form 


/[_ 

7’* 


r iC 

T* 



( 2 ) 

w here /< and 7 are the viscosity and temperature whilst a subscript w denotes a quantity 
evaluated at the wa U The basic boundary layer flow has smooth velocity temperature 
i ensi y pro i es (Stewartson [16]) and may be calculated routinely. Here we concern 
ourselves with the large amplitude waves (induced hy viscosity, which are known to occ r 
.n such flows. A I, near theory of Tollmien-Schlichting waves based on triple-deck theory 
was first given by Smith [13J for incompressible flows. Subsequently Smith [15] developed 
a inear stability for compressible flows and our concern here is with the development of 
these Tollmien-Schlichting waves in a strongly nonlinear regime. 

In order to describe the nonlinear evolution of Tollmien-Schlichting waves in a com- 

Wiflia TsT ^ ^ “I" theor ~ v based the work of Stewartson and 

Williams [18] and Stewartson [16]. We introduce the small parameter c defined bv 


c = R -'/ 8 


( 3 ) 


and look for a solution of the Navier Stokes equations for c « 1. 

Suppose that (x,g, r) and <>t. v, w) are dimensionless coordinates and velocities scaled 
on A and respect, vely. In the triple-deck limit the flow responds on a short 0(< 3 ) 
lengthscale in the X direction and the lower deck must be taken to be of order e 5 if vision] 
and convective terms are to balance. The main deck is simply the region occupied by 
he undisturbed boundary layer so that y = 0(<<) there. Finally the upper deck where 
he flow , a potential one, is defined by y = 0( f 3 ). Since the streamwise lengthscale is 

i„T hTu 7 " e T u“ enti0n ‘° S ° me "^‘‘^rhood around the position where 

variah! y v fir l°7 C | rS ' ' ,>osltl °” ™ r ''«Pand s to x = x then we define triple-deck 
variables A and Z by writing * 


x - x 


where 


- V — - 7 

./<■• ,/,S - ^ 

(4a) 



(Ml - l)3/BA«/4- 

(4b) 
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Here J is an 0(1) quantity and A/*, is the free stream Mach number and A is the un- 
perturbed skin friction at x = T. The variable Y is 0(1) in the lower deck and defined 

by (5) 


’ " f»C s /»(T;/Ti) 3 />- 
In the main and upper decks the dimensionless variables V , > are defined by 


V = 


y 




(6a) 


y y^ 5/4 ( A C ~ U' 78 (6b) 

£ 3 C 3/8( T ./^)3/2’ 

respectively. Note here that the 0(1) scaling factors have been introduced in order to 
simplify the form of the triple-deck equations which we are now in a position to write 

down. 


Lower Deck Equations 

Here it, r, i a, p - Poo expand as 

u = tC xl \r;jT^ l2 \ Xl \Ml - 1 y X/S U + 0(e 2 ), 

v = t 3 c 3/8 (r u *,/r^) 1/2 A 3/4 (M^ - ifv + o(e 4 ), 

u, = eC' 1/8 (r;/T < ;) 1/2 A , / 4 (M^-ir 1/8 fT + 0(e 2 ), 


/ = 


p- Poo = c 2 c , / 4 A‘/ 2 (AC-i)- 1/4 / J + 0(c 3 ). 

If we define a dimensionless time variable t by 

r *C(Ae>-i ) 1/4 * 3/2 

LcHWT^C'/* 

then the zeroth order approximation to the Navier Stokes equations in the lower deck 
yields 

U X + Vy + W Z = 0 , 

Ut + UUx + VUy + WUz = - Px+Vyy , (7) 

W, + VWx + VVy + WWz = -Pz + Wyy, 

Py = 0 . 

Thus the flow is determined by the 3D unsteady boundary layer equations in the lower 
deck Note also that the energy equation plays no role in the lower deck since the 
temperature and density are uniform there, see for example Stewartson and Williams 

(1969) for a discussion of this point. 

At the wall we impose the condition 


U = V = W = 0 , Y = 0 , 


( 8 ) 


whilst at the edge of the lower deck 

U-* Y + A(X,Z,T),W -*Q. 


( 9 ) 
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The linear shear flow is of course the near-wall form of the unperturbed flow and A(X, Z, T ) 
is a displacement function coupled to the pressure by the appropriate pressure displace- 
ment law obtained by solving the equations of motion in the main and upper decks. 


Main Deck 

This layer is relatively passive and transmits the displacement generated in the lower 
deck to the region above the unperturbed boundary layer. Here u,v,w and p expand as 


T* 


U o{Y ) + ~ D- 1/8 A -^AU’(Y) + 0(c 2 ), (10a) 


V = 


= — c 2 C 1/4 \'t 2 (Ml - 1 )' /4 A x Uo(Y) + 0(c 2 ), 

A oo U 0 (Y) 


(10b) 

(10c) 
(lOd) 

together with similar displacement perturbations to the unperturbed temperature and 
density. Note also that U 0 (Y ) and R 0 (Y ) are the unperturbed strearnwise velocity and 
density respectively whilst from the spanwise momentum equation the function Q satisfies 


P-Po = t 2 C l/4 * l/2 {Ml-l)- 1 ' 4 P(X,Zj)+0(e 3 ), 


Qx = —Pz 


( 11 ) 


Upper Deck 

Here the flow is potential and may be obtained by writing 

r2 £’l/4\l/2 

= [1.0.0,,,*,] + _ l ) 1 / 4 , pj + 0(€ 3 ), 

and after some simplification we obtain 

(M» ~ 1 )(pxx - Pyy) - Pzz " 0, 

P -> 0,F -» 00, 


( 12 ) 


p 


p 

' dY 


d 2 A _ 
dX 2 aS } 


0 


+ • 


The solution of (12) implies a pressure displacement law which we write symbolically in 
the form * J 

P = C(A). (13) 

Thus the nonlinear problem specified by (7), (8), (9), and (13) is now closed and of 
course must be solved numerically unless some simplification is made. Linear Tollmien- 
Schhchting waves may be described, following [15], by perturbing about the unperturbed 
flow r (J — } . Thus in the lower deck w r e write for example 


3 



and the frequency, and wavenumbers of the disturbance are found to satisfy the eigen 


relation 


m o) w p 






(14) 


where Tj = — ? afto T /b is the Airy function and 

Ai(q)(lq. 

The eigenrelation (14) can be solved for a complex wavenumber a for given values 
of /?, M '*> and ft. The flow is unstable or stable depending on whether a, is negative or 
positive" If we fix f) and AT* we find that -a, increases from a negative value when 
n increases from 0. The growth rate then changes sign at ftyv, the neutral value of 
ST and then increases monotonically until ft = ft A f where the maximum growth rate 
is" achieved. Beyond ft = ft M the growth rate can have further maxima of -a, but 
ultimately asymptotes to 0 at large frequencies. Note however that the most unstable 
wave always corresponds to the first maxima when & = 0 and that -o, remains positive 
for ft >> 1. We can alternatively choose to keep a real and find the complex value of U. 

For weakly nonlinear Tollmien-Schlichting waves it has been shown by Smith [14] and 
Hall and Smith [7] that nonlinear effects are stabilizing for both two and three dimensional 
waves. For axisymmetric compressible flows Duck and Hall [3], [4] showed that nonlinear 
effects are again stabilizing but we know of no results for planar compressible flows. The 
analysis we shall give is appropriate to strongly nonlinear disturbances but in the small 
amplitude limit of our theory we will be able to infer the role of nonlinear effects in 
the weakly nonlinear regime. Strongly nonlinear waves are governed by the full triple 
deck problem (7)-(8),(1.9), and (13). In general such disturbances must be described by 
numerical means, here we shall investigate asymptotically large amplitude locally periodic 
wavesy stems. 



2 The phase equation for large amplitude waves 

The essential ideas of the approach we use can be found in Whitham [21] and were pre- 
viously used by Hall [6] to discuss the two-dimensional incompressible case. We suppose 
that a wavesystem moves through the flow and that the local streamwise wavenumber and 
the frequency of the waves vary slowly as the wavesystem moves downstream. Suppose 
that the slow temporal and spatial variations of the wave occur over times and lengths 
of relative size compared to the period and wavelength of the wave. We introduce 
slow variable and r by writing 

£ = 6X, (15a) 

t = H, (15b) 

where 0 <<!><< 1. In order to represent a locally periodic wavesystem we define a phase 
function 0 by writing 

e = -Q(^T) + liZ. (1G) 

V 

Here ft is a constant spanwise wavenumber and the streamwise wavenumber a and 
the frequency ft are then defined by 

O = 0 O ft = -0r, 
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( 18 ) 


so that, n and Q must satisfy the evolution equation 


o> + n* = 0. 

The frequency !! = !!(«,«) may then he found in terms of an asymptotic expansion 
involving 6 and in principle tile evolution equation ( 18) can be solved to determine how a 

vanes ,n tune and space. The definition (2.1)-(2.3) imply that in the triple-deck equations 
we must, use the transformations. 


() 

dX 


« 


<L , cA l 

<90 dt'dt 


n d <) d 

0— -f 8 — , — 
d0 dr dZ 



The lower deck solution for y = 0(1) 

We expand U, V , and W in the lower deck in the form 


(t/, V, H ) — (Vo, Vo, H o) + (U\, \\ , H | ) + . . . ( 19 ) 

where Vo,V 0 , etc. are functions of 0,X,V and r. In (19) we have anticipated that the 
appropriate choice of expansion parameter is 8J rather than 8 ; the reason for this choice 
will become clear when we investigate the form of the leading order terms in (2.5) in the 
limit y -k oo. The pressure is then expanded in the form 


P = 6 2/3 A/(£, T) + P 0 + 8' /3 f\ + . . . 


( 20 ) 


where P P, etc. are functions of $,(, V and r whilst P„ is an induced pressure again 
implied by the form ol ((/„, V„, H oi for ) » 1. The frequency then expands as 

=: Do + 8d 3 Qj 

The leading order approximation to (1.7) is then found to be 


otl'oe + Vov + fiWo s = 0 

-MoVoe + al> 0 f l oe + V 0 Uo Y + dW 0 Vo$ = -aP oe -f U 0YY , 

-tooWo* + aUoWoe + V 0 H’ 0 > + PWqWw = ~f3P o0 -f W 0YY , ( 21 ) 

/ J ov = 0, 


and we seek a solution of these equations which is periodic in 0. 
function A is expanded as 

A = A 0 + £ 1/3 d, + . . . 

then the conditions required to close the leading order problem are 

Uo = Vo = Ho = 0, V = 0, 

Vo — * V + c4 0 , V — > oo, 

Po = C(Ao). 

The partial differential system (2.7)-(2.9) specifies a nonlinear eigenvalue problem for fi 0 
as a function of a, ft and M v . We note that the disturbance size can then be changed 
as a function of o and that solutions will only exist for certain values of a. Moreover it 
follows that the slow dependence of (f/ 0 , V 0t H o) and P 0 on r arises only through the 


If the displacement 


( 22 ) 
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dependences of these functions on o and fi 0 . In order to solve (2.7)-(2.9) numerically it 
is convenient to write 

oo 

(Co, Vo, Wo, Po) = (Coo, 0, Woo. Poo) + Y (Co.<, W 0n , Won, Po,n)P (24) 

n^O 

After a little algebra it is possible to show that the partial differential equation (2.7) can 
be reduced to 

“ in(«r + nf'oo + - SloJ'o" + >Von{aCoo + ^W-oo}" = 

V -2— {C-mVihn-Vi',^0 (25) 

n — m 

m^n 

m^O 

for n — ±1, ± 2 , . . ., 

o l qq + /? Wqo — Y* ~~ { Wom V 0ni — W om W 0m } , (26) 

m# 0 m 


Coo — Woo — Vc 


On 


— V" 

— V 07l 


= 0 , n — ± 1 , ± 2 , . . . , 


' On 


ffi — w 1 =0 
1 oo — * K oo u ’ 


Y = 0 , 
V* — * oo . 


(27) 


C( 0 ) = 


mo (a 2 + /? 2 )V'o„(oo) 


\Ml - 1 Q ) 


1/2 


(28) 


The last equation is obtained by projecting (1.13) into Fourier space and relating V^(oo), P 0n 
to the corresponding displacement function Fourier component (through continuity) and 
Co„, Wo' n at the wall respectively. Of particular importance is the fact that the mean flow 
correction terms only appear in the leading order nonlinear eigenrelation in the combi- 
nation oCoo + /WVoo so that we do not need to solve for Uoo and Woo independently. In 
addition we see that the mean flow functions tend to constants at the edge of the lower 
deck. We write 

Coo - * Aw? Woo —+ Poo? V * 00 1 


and note that /loo, Poo must in fact depend on the slow variables £ and r through the 
dependence of these quantities on a. This is a crucial result because it means that the 
mean flow correction has an outer layer structure where ~ 977 • The outer layer 

is therefore of depth 3 and it is for that reason that our lower deck expansions were 
previously taken to be in powers of 6 l/3 rather than S. We point out at this stage that the 
functions Coo, Won, W 0n , Po„,» # 0 do not develop such an outer layer structure because 
they possess a fast dependence on A through the phase function 0. 

Before we discuss the solution in the outer < 5 1/3 boundary layer it is convenient to 
write down the order approximation in the Y = 0 ( 1 ) part of the lower deck. The 
appropriate differential equations are found to be 
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(30a) 


! u> + V\Y 4 - ftW\8 — 0 , 
-n 0 r ie + a{H o U l0 + IhVoe] + Vo UlY + \\Uoy + 
PiWoVio + Wd'oe) + ot P\ e - V XYY = tt, V 0 o, 
-n 0 Wie + a{U 0 W ie + UiWob} + V X W 0Y + V 0 W tv + 
/4{H 0 H 4- H 1 1 1 (w } -f [3 I\ff — H'n y = H 0(? , 

PlY = 0. 


(30b) 

(30c) 

(30d) 


The m a,n pent to not.ee here ,s that {U u V U W U P X ) = &«/ 0 , Vo , W 0 ,P 0 ) is a solution of 
the homogeneous form of these equations. This is because of the invariance of the 0(1 ) 
solution to a Ga hlean transformation in the A direction. The inhomogeneous form of 
(z.16) is to be solved subject to 


Cl = Vi = IT, = o, V = 0, 


Px 


(31a) 


(31b) 

and a condition on U X (V = oo) obtained by a consideration of the outer boundary layer 
The condition that the resulting partial differential system should have a solution will 
determine SI, as a function of o. This enables us to write down the approximation to the 
phase equation (18) correct upto order 6 1 ' 3 . I n the absence of the outer layer we would 
wnte < °" n the c °nditions C, - A u W x -> 0, Y -> oo. However the required conditions 
can only be written down once the solution in the 0(6‘/ 3 ) diffusion layer has been found 


The diffusion layer 


The dominant balance in this layer is between diffusion in the V direction and convection 
in the A direction. The thickness of the layer is governed by the balance U& ~ -2L 

V S n',T./ 3 , ~ ’ fOT ''>> '• a,ld W ~ m <» follow, that the required llyer hU 
y ~ 0(6 */ J ) so we introduce tj defined by 

i) = 6 l '*y. 

We now expand C, V, W and P in the form 


V = {6 l/3 i) + U M {t ]) + 0(<5 1/3 } + [f/ 0 + S l/3 U X + 0(6 2/3 )j 

V = {6 2/3 Vm(ii) 4- 0(6)} 4- V 0 4- 6 1/3 \\ + 0(6 2/3 )J 3 , 

W = { W M (, f ) 4 - 0 ( 6 1/3 )} + [l| 0 + 6 ,/3 li / , + 0 ( 6 2/3 )J a'/ 3 , 

P = {6~ 1/3 p M + 0(1)+} + [p 0 + S 1/3 P X + 0(6 2/3 )J . 


(32a) 
(32b) 
(32c) 
(32d) 

In the above expansions the terms in the curly and square brackets correspond to mean 
flow correction terms and the wavelike part of the disturbance flow respectively If we 
substitute the above expansions into the rescaled form of the lower deck equations and 
match with the solution in the } =0(1) region it is easy to show that 

I 0 = Aq — A 00 , Vo = 
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The mean flow functions V M , V m , W m are determined by the solution of 

ijl \i( T Vm = t I'Mipii 

Vm t + ^A/rj = Oi 
Vm = ^ooi 7 1 = 

Vm ~~ * Am{£,)iV 00 ’ 

Pm = ^(Aa/), 


(33) 


and 


= W'Mtjt,, 

H a/ = Boo-, 7 ! = 0’ * 

11\/ — + 0, 7/ — ► OO. , 


(34) 


These systems are most easily solved by taking a Fourier transform with respect to £ 
and indeed (33) has been solved previously in the context of linearized triple deck theory 
by Stewartson [17] and Smith [12]. For our purposes here we only need the quantities 
( T/jp IV'a/t) evaluated at 7 / = 0. We obtain 


V' M (0) = Mi,?) = 


_3A l (0) S{i _ 

Ztt J —00 


VV r ; f (0) = J 2 (e,r) = 


/T 1 1 

V 277 3 , / 3 A,(0) 



1 

(Z-<t>yv 




(35) 


where 



V’ 4/3 (l + V’ 1/3 )e~ flv ’ 

1 -f j/> 4 / 3 + j/’ 8 / 3 


dV\ <j> > 0, 


S(d>) = 



,d> < o 


with 

0 = — 3i4'(0). 

The functions f'y, Vi and IT, are then found in terms of V M and W M . We obtain 

Vj = —qAo$Vm — fiAo$W m — a Qo(i ( 36 ) 

aVi + fiW\ = {A 0 - Aoo)(aV M + pW’ M ) + Qoa. (37) 

where Q 0 is to be determined by matching with the main deck solutions. When 7/ — ► 00, 
V' m ,W'm — 0 so that Q 0 = At. 

It remains for us to match (2.22), (2.23) with the solution in the Y = 0(1) part of the 
lower deck. In the first instance we note that Vm ~ JiV,Wm ~ J 2 V f° r sma ^ 0 so ^ ' s 
convenient to write 

V\ = V\ + J\ 7/, Wi = W , 4- J 2 7/, 


in (2.16) to give 


aV\o + V 1 v + /3VV'i e — 0? 

-floVio + q(VqV\9 + (hV 0 g) + V 0 V\y + V X V W + Pi WoVie + W\Voo) 
+aPie - V\YY = 00 - ot-hyVoe - Vo J, - P-h i)Vm 

-HoVVie + a{U 0 Wie + ViWoo) + V\V oY + V 0 \ V lY 
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(38) 


+ 0{W o W, e + W x W m ) + (iP ie - Wiry 

= ih Woe - m/./, W Q e - V 0 J 2 - flJwWoB. 

The conditions required to completely specify the problem for (U t ,V\, W u P t ): are then 
found to be 

l\ = W\ = V x = 0, V = 0, 


Pi=C{A 1 ,), 

- qA\ + aUj + fiWi — » (A 0 - /l 0 o)(«-./i 4- /?J 2 ), > — > oo. (39) 

Suppose that we set all the right hand sides in (2.24) equal to zero and replace the above 

condition by o( r , + fi\\ x -» qA u ) -> oo. The simplified homogeneous system obtained 
by this procedure has the solution 


(Ui,V u W u Px) = (U 0 ,V 0y W 0 ,P 0 ) e 


so that (38), (39) only has a solution if the appropriate solvability condition is satisfied. 
The condition is most easily obtained by expressing (38) as a system of equations in the 
form 

d 


d\ 


-■(\.V,. \v./V) 7 = B{x,V u \y,P x ) t + C(x,Vuxy,P\)1 + G, 


where, \ = aV x + fSW x , B,C are 4x4 matrices and G is a column vector with elements 
proportional to aJ Y + f3J 2 . The partial differential system adjoint to the homogeneous 
form of the above system is obtained by multiplying by a vector and integrating over 
0 < 4 < oo, 0 < 0 < Full details are given in [6] and after some analysis we obtain 


fii = I\'{a)(Qj l + /iJ-2 ) (40) 

where A (a) may be expressed in terms of integrals involving the adjoint and the 0(1) 
lower deck solution. 

Having determined flj we can now write down the partial differential equation to 
determine o correct to 0(<S'/ 3 ). We write 


■ )t + + ^ 1/3 {A'(o)(qJ, + = 0. (24) 

If o is then specified at r = 0 we can integrate (24) forward in time to determine the 
evolution of a. 


3 The modulational instability of uniform wavetrains 

It can be seen from (24) that a particular solution of the wavenumber evolution equation 

IS 

a = o 0 

where o 0 is a constant so that the wavenumber and frequency of the wave remain constant 
as it. propagates downstream. In order to investigate the modulational instability of the 
uniform wavetrain we write 

a = o 0 + A 


9 



with |f-| << 

’ Of 0 ' 


1 . It, follows from (24) that A satisfies 


~ + [%{ Qo ) + MV'(ao) + 0( A 2 )] ^ + 6'/ 3 K(a o)^ [aJi + M] - 0(<$ 1/3 A 2 ). 


Here J\ and J 2 are defined by 


Jo 


1 Ko(*°) j* dA 


1 


-M. 


2 tt 3 1 / 3 4,(0) J-oo^U -<^) l/3 

Bv a suitable transformation we can eliminate the term in the wavenumber evolution 


DV ct suitauit - • t 1 c 1 /3_ 

equation proportional to Q ' 0 (a 0 ). If we then introduce a scaled tune variable > ~ 4 r ’ 
and amplitude A ~ 0 (r 1/3 )A then we see that the canonical form for the evolution 

equation is ^ r* d \ { H (£ - <f>) „ 

ar + A ^ = ± ^i- 

Here H is the Heaviside function and 7 is defined by 


■1 


4>) f d(f> 


(25) 


7 


W^_L. 

Oot4nn(oo) 3 l/,<5 


(26) 


Calculations show that §«• is always negative so that 7 is a positive constant. Finally 
the ± sign in (25) are to°°be taken depending on whether B' 00 (o 0 )K{ao) is positive or 
negative. Computations which we shall report on in the next section show that both 

possibilities occur. . ... 

In the absence of the viscous like terms on the right hand side of (2o) we see that 

A satisfies the inviscid Burgers equation which for rather arbitrary initial data will de- 
velop shocklike solutions. The viscous terms will then become important in such rapidly 
varying regions. The kernel function S appearing in (25) can then be simplified using its 
asymptotic form for small argument. Thus using the result 

S{<{>) ~ ^-T( l/3)<r 1/3 +,0 < <f> << 1. 


s 


37T 


2\/2tt 

we deduce that a rescaled form of (25) is 


,<i> < o, 


dA 

OT 


+ A 


dA _ t< 
df ~ dU- 


dZ 


d_ 

di 


r( dA 


1 


d<t> (£ - <£) 1/3 


de 


and we shall now concentrate our attention on this equation. 


(27) 
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4 Numerical solution of the 0(1) problem 

Here we discuss the results we have obtained for the solution of the nonlinear eigenvalue 
ptoblein Q 0 — D 0 (cv, /3, A/ x ,) associated with the 3D periodic triple deck problem (2.7)- 
(2.9). I he calculations were carried out using a finite difference method in the Y direction 
and a houner spectral method in the 9 direction. The nonlinear terms were accounted for 
by an iterative process with only the linear terms evaluated at the new level of iteration. 
We restricted our attention to three Mach numbers, A/ x , = 1.1, 1.25, 1.5. Figure 1 shows 
the dependence of o on ft for the linear neutral modes at these Mach numbers. The 
corresponding frequencies are shown in Figure 2. 

Figures 3-5 show a variety of Q 0 - o plots for three different values of A/ x , and a 
selection of spanwise wavenumbers. The circle shown in the figures corresponds to the 
linear neutral point in D 0 - «o space. Thus one end of each of the curves connects the 
zero amplitude state with an initially small amplitude wave whilst the other end of each 
curve represents the points at which our iteration procedure failed to converge. In each 
case we used 32 modes in the periodic variable 9, 400 in the normal direction and infinity 
was taken to be 15. At the point on each curve where the iteration procedure failed to 
converge we calculated the wall shear stress but there was no sign of it changing sign. 
We have no explanation why our iteration procedure could not be continued further and 
it is possible that, as was found in [6] for the incompressible case, the Fourier expansion 
procedure fails to converge at this stage and the disturbance develops a singularity. 

Figures 3a-i show results at a Mach number 1.1. At moderate wavenumbers, Figures 
3e-g, we see that there is a supercritical bifurcation to a finite amplitude state at the 
linear neutral point. 1 his is similar to the situation discussed for subsonic flow bv Hall 
( 199o). However in f igures 3b, h we see that in the weakly nonlinear limit the wavenumber 
initially decreases as the amplitude increases. Since the frequency initially increases 
in this regime the bifurcation remains a supercritical one with the bifurcating solution 
temporally stable. However for sufficiently small or large wavenumbers. Figures 3a,i 
we find that both the frequency and wavenumber decrease with increasing disturbance 
amplitude in the weakly nonlinear limit. Thus the bifurcation is now subcritical and 
nonlinear effects are destabilizing. 

Thus all disturbances apart from those with sufficiently small or large spanwise 
wavenumber bifurcate supercritically. The modes with moderate spanwise wavenumber 
bifurcate subcritically. 

In Figures 4a-f we show results for a Mach number 1.25. The results are similar 
to those described previously for A/ x = 1.]. However we see that the supercritical 
ehaviour is now restricted to a smaller band of spanwise wavenumbers. If the Mach 

number is further increased we ultimately reach a stage where the bifurcation is always 
subcritical; see Figures 5a-e. 

finally it remains here for us to discuss the sign of the constant fi' 0 (o 0 ) A'(o 0 ) which 
m effect determines whether the viscous-like derivative appearing on the right hand side 
of 3.1) is stabilizing or destabilizing. We recall that in in the subsonic case it was found 
in [6] that viscous effects were always found to be destabilizing so that the negative sign 
in (3.1) was to be taken. In the present problem we found a similar result for every 
case when the initial bifurcation has both frequency and wavenumber increasing with 
the disturbance amplitude. In every other case we computed we found that viscous 
effects were stabilizing if the equilibrium value of a initially decreased with disturbance 
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amplitude. Thus it is necessary for us to discuss the solutions of (3.1), 01 in fact more 
appropriately (3.3), for both cases. 

5 Numerical solution of the phase equations 

In this Section we are concerned with the numerical solution of the canonical phase 
equation (3.3). The integral term on the right hand side is a iractional derivative of order 
4 aiK | is most easily dealt with in Fourier space. A fractional derivative of order p is 

defined by 

l—f R. - dd>. 

dp ~ r(i -p) J - oo (Z-w 

For a discussion of fractional calculus the reader is referred to [10]. The Fourier transform 

giVeS (9 A\ 

^(^)=(^ (A) , 

where jF(u) = ue~ iki d(. According to a rescaling of A similar to that of Section 3, 
the canonical equation (3.3) is cast into the form 

d^u 

u t + mix = ±"^7 3’ 

(,r,0eR 1 xR + , (* 28 ) 

u(;r,0) = u 0 (;r),w(.r + 27r,f) = u{x,t). 

Here „ > 0 and can be re-scaled to unity by changing the initial condition, for instance. 
The boundary conditions in (5.1) indicate that we are considering the spatially periodic 
problem since our main interest lies in the prediction of local structures such as finite-time 
infinite slope singularities, for example. Equations containing fractional derivatives have 
been derived and studied by other authors in the context of boundary layer stability, [1], 
and nonlinear acoustic waves, [19]. The former work derives an equation with negative 
diffusion and a ± spatial derivative, while the latter study considers a positive diffusion 
and a fractional derivative of ± together with the usual stabilizing Burgers diffusion 
(second spatial derivative and positive diffusion). The former authors conjecture and 
formally analyze the formation of a “shock” singularity driven by the inviscid Burgers 
equation, to leading order. Our numerical results for the present problem lend support 
to such breakups. The analysis and numerical experiments in [19] do not produce shocks 
due to the presence of the stabilizing second derivative term; this is brought into the 
evolution as a higher order correction once a shock begins to form. In the present problem 
a hierarchy of higher derivatives enter simultaneously and consequently the evolution in 
the “inner” region is governed by the full unsteady triple deck problem which is hard to 
analyze in order to obtain jump conditions, for instance. 

A linearization of (5.1) leads to normal mode solutions 

u = Ae ±^)' l3 t + c.c = Ae**' + c . c ., (29) 

where c.c. denotes complex conjugate. We see from (5.2), therefore, that the system is 
linearly stable or unstable respectively, depending on the sign chosen. The linearly stable 
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case provides a nonlinear evolution problem which can be compared with the viscous 
Burgers equation. There are two essential differences, however: First the amount of dif- 
fusion is less (for large k it is of order k‘ i/:i as opposed to k 2 for the Burgers equation), and 
second, the linear solution contains a dispersive part in contrast to the purely diffusive 
nature of the linearized viscous Burgers equation. The latter fact requires some care in 
the numerical treatment described later. The main question we address numerically is 
whether the dissipation in the present problem is sufficient to prevent infinite slope singu- 
larities - the numerical evidence suggests that shocks do tend to form. The situation with 
the negative sign in (5.1 ) is worse in that the problem is now ill-posed and solutions may 
exist for small enough times if the initial conditions are chosen appropriately. This class 

of numerical experiment also suggests the generic formation of infinite slope singularities 
or “shocks”. 

The numerical scheme used is a pseudospectral one which allows a straightforward 
representation of the fractional derivative. All derivatives are computed by a forward 
Fast Fourier Transform (FFT) along with the appropriate inverse FFT, and the time- 
marchmg is clone in real space by an adaptive second or third order accurate scheme. In 

order to monitor the accuracy of the solution as the computation evolves, we utilized the 
conserved quantity 

=0. (30) 

Most numerical experiments have initial condition u 0 (x) = sin(a-) and so / 0 2ir u(.v,t)dx = 0 
throughout the evolution. All results reported here satisfy (5.3) to within computer 
round-off errors. Before presenting numerical results we mention some accuracy require- 
ments of the numerical schemes. In view of the linear solution (5.2), it is seen that when 
the *th Fourier component is stepped forward in time by an amount At (according to a 
tune-split, scheme for instance - this is a useful simplification to get an estimate of the 
accuracy requirements) the solution due to the linear operator is given by 

Uk(t + At) = ,=Mi/2)(-i+.V5)*«/ 3a l{ik(t) (31) 

If a computation is performed with n modes, then the maximum available Fourier compo- 
nent is k mnr = f (in practice k maT is smaller than n/2 due to filtering); the complex part in 
the exponential in (5.3) causes a rotation which will be spurious unless ft = I ^k^ At 
.s small enough. This places a restriction on the value of At. All runs reported here have 

/y < _° u 1 ’ a,ld 11 IS seen t!mt the ste P size Iiep ds to be decreased as the resolution increases 
lhe first set of results presented below have u = 0.1 and an initial condition t/ 0 (.r) = 
sin(.r). All numerical experiments terminate in infinite slope singularities after a finite 
time and so a “shock” is seen to form then. The maximum absolute value of the slope 
that can be computed accurately depends on the spatial resolution and it was found 
that an amplification in ma.r(|u T |) of over 80 is possible to achieve accurately with 4096 
modes. Numerical convergence has also been established by a comparison between this 
run and a lower resolution run. Unless otherwise stated the results shown in the following 
figures have been computed with n = 4096. 

Figure 6(a) shows the evolution of «(*,/) ami Figure 6(h) depicts the final computed 
profile. 1 he maximum value of the magnitude of the slope is 80.03 by the end of the 
computation, and it occurs at t = 1.246. It is seen that the solution steepens up as the 
computation evolves and due to the dispersive part of the pseudo-differential operator 
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the position of maximum |u,| is not constant but varies slightly. Next we consider the 
evolution of the maximum magnitude of u T in Figure 7 which again indicates that the 
slope is blowing up after a finite time. We note that the initial value of this quantity is 
1 and by the end of the computation a magnification factor of 100 has taken place. The 
energy of the solution defined as jf u 2 (xj)dx remains bounded as t increases and in fact 

decreases monotonically; this is shown in Figure 8. 

Further quantitative information about the suggested singularity can be obtained by 
considering the evolution in time of the spectral amplitude of the solution. Given a 
time t , say, the discrete Fourier transform of u is computed and is denoted by u* for 

f, _ ii, . . . , ” ; it is sufficient to consider those spectral components with positive k 

(note that the component corresponding to k = 0 is zero due to the conservation of the 
integral of u.) since u(;r,<) is a real function and so the sequence {u k } is Hermitian, i.e. 
u_ k = vl where * denotes complex conjugate. We define a spectral amplitude, therefore, 

at a given time t by n 

S(k,t) = \u k (t)l k = 1, — 1- 

In Figures 9 and 10 we depict S(k,t) versus k and In {S{k,t)) versus In (k) at different 
times. The times correspond exactly to the times at which the profiles in Figure 5 are 
shown and are therefore the spectral amplitude characteristics of those profiles. In both 
Figures 9 and 10 as time increases the spectral amplitude increases upwards. In order 
to capture the behavior of the spectrum for large k we follow [20] and assume that the 
spectrum has the following behavior for large k , 

u k (t) = c(t)k- o{t) e- i{l)k . ( 32 ) 

The form (5.5) is suggested by the large k behavior of the Fourier transform of an ana- 
lytic function with isolated singularities in the complex plane. The dominant behavior 
contained in (5.5) comes from the singularity in the upper half-plane nearest to the real 
axis which is not a multiple pole (for details see [20] and [2]). We used the results of 
Figures 9 and 10 to estimate the parameters in (5.5) as the computation evolved towards 
a singularity. We note that 8(t) decreases to zero as the singularity is approached and the 
computation is inaccurate once the value of 6 becomes smaller than the mesh size (this 
accuracy requirement is observed by all the results given here). We used least squares 
fits to estimate a(t) and 8(t). It. follows from (5.5) that 

In |ttfc| = — a{t) \n(k) - k8(t ) + I)(t), (33) 

where, for notational convenience, o and 8 in (5.6) are real numbers which correspond 
to their respective real parts in (5.5). It follows from (5.6) that log-linear and log-log 
plots will provide estimates for 8(t) and a(i) respectively. The former evolution can be 
found in Figure 9 and the latter in Figure 10. The results for our estimates are given 
in Figure 1 1 and summarized in Tables 1 and 2 which contain additional information on 
the modes retained for the fit. The spacing between grid points for 4096 modes is 0.0015 
and is larger than the minimum value of 8 reached by the end of the computation. The 
numerical results can be used to conjecture the terminal form of the solution. It appears 
that the dominant structure is that provided by the inviscid Burgers breakdown - the 
inviscid Burgers has a power law dependence of a = 4/3 due to the local :r behavior 
in the neighborhood of the singularity. At times above t = 1.24 a power law of 1.33 is 
achieved suggesting that it is the inviscid Burgers mechanism which is controlling the 

local structure of the singularity. 
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t 

a 

Modes retained for fit 

0.7 

3.03 

2 < k < 8 

.8 

2.61 

2 < k < 10 

.9 

2.38 

2 < k < 20 

1.0 

2.06 

5 < fc < 30 

1.1 

1.76 

10 < k < 60 

1.2 

1.42 

10 < k < 100 

1.21 

1.39 

10 < k < 100 

i .22 

1.37 

10 < k < 100 1 

1.2.3 

1.35 

10 < /• < 100 

1.24 

1.33 

10 < k < 100 

1.246 

1.33 

10 < k < 100 


Table 1 


t 

6 

Modes retained for fit 

0.7 

0.377 

30 < k < 70 

.8 

0.168 

40 < k < 150 

.9 

0.099 

50 < k < 250 

1.0 

0.051 

50 < k < 490 

1.1 

0.021 

500 < k < 1150 

1.2 

0.0082 

1000 <k< 2000 

1.21 

0.0076 

f 000 < h < 2000 

1.22 

0.0071 

800 < k < 1800 

1.23 

0.0067 

800 < k < 1800 

1.24 

0.0063 

800 < k < 1800 

1 .246 

0.0060 

800 < k < 1800 


Table 2 

Another quantity which provides positive evidence of an inviscid Burgers breakdown 
is given by the enstrophy of the solution defined by 

2£(t) = f u 2 x dx. (. 34 ) 

The evolution of £{t) corresponding to the runs described above is given in Figure 12 . In 
order to make a comparison with the inviscid Burgers analogue, it is useful to give the 
solution then, corresponding to i/ = 0. Using characteristic coordinates an implicit form 
of the solution is (see [ 9 ]) 

w(i,ar(f,s)) = «o(s), 
x(t,s) = s + tu 0 (s). 

The slope u r follows from these results and is 


1 + t u'q( s ) 
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showing that an infinite slope singularity is encountered after a finite time given by 

1 

T ~ -Inf[u' 0 (s)}' 


An estimate of the rate at which £(t) blows up as the singular time is approached is best 
achieved by using Lagrangian coordinates; this is described in [5] and we simply sketch 
the result. In Lagrangian coordinates the enstrophy is 



K(s)) 

1 + tU()(s 


If the singular time is <„ then it follows from a local expansion of u 0 near s - 0 (without 
loss of generality the singularity is taken to occur here) that 


which on insertion into (5.7) and retainment of the leading order term as t -> L , gives 


£{t) ~ (tm — 0 2 ■> t * 


(35) 


The local r 1/3 structure of u near the singular time also follows from the Lagrangian 
formulation (see [5] for more details). The results given above are valid for the inviscid 
Burgers equation and in general the singular time U is different from that of (o. 1 ) denoted 


In what follows, then, we show that our numerical results with v # 0 are consistent 
with the asymptotic blowup rate (5.8). Referring to Figure 12 we see that £ is becoming 
singular after a finite time; the value of t. is not provided directly by the computations, 
and in order to check the algebraic rate (5.8) it is necessary know the value of t s . In the 
absence of a rigorous wav of obtaining t s we proceed as follows: reasonable guesses are 
made for consistent with the data in Figure 12, and plots of ln(£(t)) versus ln(f , - t) 
are generated (for brevity we do not include those plots here). If a power law behavior is 
present for small values of (t, - <), then the plots above are linear with the slope giving 
an estimate of the power. Least squares fits were applied for a range 1.226 < t < L24b ( a 
typical corresponding range in the log-log plots is [-2.9, -2.59]) to gain an estimate of the 
power law for comparison with (5.8) which has a power law behavior of — The resu ts 
are given in Table 3 below. The numerical values are consistent with those predicted by 
the local analysis providing additional evidence that it is the singularity formation of the 
inviscid Burgers equation which is dictating the singularity structure. 


t s 

1.29 

1.3 

1.301 

1.3025 

1 .305 

1.31 

slope 

-0.4P 

-0.49 

-0.50 ' 

-0.51 

-0.53 

-0.57 


Table 3 


5.1 Negative diffusion case 

In this section we consider the evolution when the “diffusion” is negative and the minus 
sign is in effect in equation (5.1). In view of the linearized solution (5.4) the problem 


16 




now contains a short wave instability with growth of the order exp (k 4 ' 3 th we note that 
this is worse than the well-known Kelvin-Helmholtz instability whose nonlinear evolu- 
tion has been treated numerically by [8], for example. Due to the short wave instability 
numerical round-off errors are viewed by the system as low amplitude high wavenumber 
perturbations; the shorter the wavelength the faster the initial (exponential) growth will 
>e and unless such perturbations are removed their eventual growth and nonlinear inter- 
actions contaminate the flow and produce a spurious numerical solution. This problem 
is overcome numerically by filtering out (setting to zero) any Fourier coefficient below a 
certain olerance level; we use the same filter for all wavenumbers in much the same way 
as was done in [8] for Kelvin-Helmholtz instability and in [11] for a system with a worse 
dispersion relation than Kelvin-Helmholtz flow. This filtering technique was also applied 
to the positive diffusion case described earlier. 

The presence of the short-wave instability, then, means that nonlinear solutions will 
exist for small enough times if the initial conditions are chosen appropriately. To illustrate 
this consider the evolution equation with the bad sign 


Uf -|- U Uj 


— v 


dx 4 / 3 ' 


v > 0 . 


(36) 


oven an initial condition at t = t 0 < say, a backwards integration in time is well-posed, 
to t.mes prior to any singularity encounter at least (care is required since the positive 
diffusion case can encounter an infinite slope singularity if the solution is evolved long 
enough), since the diffusive effect appears positive then. Initial conditions which satisfy 
(oJ) can thus be generated. Typically we integrated (5.9) from an initial condition «„(*) 
at /„ to t - 0 and there switched to a forward time integration. In the numerical results 

cescimec )e ow, we have t„ = 0.05 and u 0 = sin(.r). The number of modes used is 4096 
arid the filter level is set at 10~ 2 . 

The evotai 0 " Iron, / = 0 to i = 0.91 is give,, i„ Figure 13. Profiles are plotted 
eve,j 0.05 tm,e units and the profile at ( = 0.91 is include,! also. Again we see that 
the solution is steepening and an infinite slope singularity appears after a finite time- 
the minimum slope at t = 0.91 is approximately -20. In addition, the maximum (or 
minimum) of the solution increases with time as opposed to the behavior for positive 
n described previous y. In fact the energy of the solution increases monotonically now, 
as expected. These features of the evolution are given in Figure 14 which depicts the 
evolution with time of the maximum of |u r |, the enstrophy (defined earlier) and the 
energy. These results strongly suggest that the numerical solution is encountering an 
infinite slope singularity after a finite time. We have also checked the power law (5.8) for 
n e enstrop y with favorable results; it is found that a singular time between 0.93 and 
■ gives the power law (5.8) for our numerical data of Figure 13. We also note that 
the numerical solution becomes unstable soon after t = 0.91; in particular we were not 

able to achieve with desirable accuracy, the minimum slopes of -80 or less which were 
achieved for positive diffusion. 

Finally we consider the evolution of the spectra of the solution. As before, we present 
log-linear and log-log plots analogous to those of Figures 9 and 10. Results are plotted at 
time intervals of 0.05 for 0 < t < 0.65 and at intervals of 0.005 for 0 65 </< 0 91 The 
reason for this refinement is due to the rapid change of the spectrum as the singularity 

is approached. I lie log-1, near plots are given in Figure 15 and the corresponding log-log 
ones in Tigure 16. 666 
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As can be seen from Figures 15 and 16 there are some different features in the spec- 
trum from that corresponding to positive diffusion. In particular an oscillatory behaviot 
appears at large wavenumbers and a kink in the spectrum is clearly seen at a time just 
below 0.91. Such features were absent in the positive diffusion case and can be explained 
by the influence of higher singularities moving down in the complex plane towards the 
real axis; the interaction of several such singularities can produce oscillatory behavior in 
the spectrum at large k. For an explicit example of this the reader is referred to the text 
[2] (p. 256). We also tried to estimate the power law behavior of the spectrum in order 
to compare with the positive diffusion results. Such an estimate can be derived from the 
data, of Figure 16 by the same techniques described previously. It has been found that for 
the times up to which the integration has been carried out, the best estimate of the slope 
coming from intermediate values of \n(k) in Figure 16, is approximately -1.78; this is 
approximately 25% different from the estimate found earlier and for the inviscid Burgers 
behavior. The solution suffers from short wave instabilities immediately after t — 0.91 
(these can be seen creeping in at the largest values of k already) and it is possible that a 
/t-dependent filter can increase the time of integration and allow further convergence in 
the power law estimate. We emphasize, however, that the estimate -1.78 is a consistent 
transient with the inviscid Burgers singularity and further careful numerical experiments 
are required to verify this. 

6 Conclusion 

We have investigated in detail the solution of the modulation equation for large amplitude 
waves in a compressible boundary layer. The equation in question, (3.3), was solved foi 
the cases of both negative and positive diffusion. The evolution equation describes the 
evolution of a wavenumber perturbation to a uniform wavetrain. In both cases we found 
that with periodic initial data the solution breaks down and a singularity forms. After 
the onset of the singularity the full triple deck problem must be solved. However our 
results show that large amplitude waves in supersonic flows cannot persist with constant 
wavenumber and frequency. 
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Figure 1 The dependence of the neutral linear wavenumber a on the 
wavenumber 0 for AT* = 1.1, 1.25, 1.5. 
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Figure 3c, beta=0.4 
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Figures 3a-i The dependence of JIq on a for the nonlinear states with = 
and ft = .1, .2, .4, .6, .8, 1., 1.2, 2.6, 3.0. 
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Figure 4a, beta=0.4 
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Figure 4b, beta=0.6 
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Figure 4c, beta=0.8 
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• Figures 4a-f The dependence of on o for the nonlinear states with M 
and B = .4, .6, .8,1., 1.2, 1.8. 
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Figure 5a, beta=0.4 
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• Figures 5a-e The dependence of on a for the nonlinear states with M, 
and (3 = .4, . 6 , .8,1., 1.2. 







Evolution of u, n=4096, nu=0.1 , t=0->1 .2 
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• Figure 6 (a) Run 1. Evolution of u(x,t ): Positive diffusion case is = 0.1, initial 
condition u 0 (x) = sin(x), 4096 modes used; profiles shown every 0.1 time units, (b) 
Final computed profile at t = 1.246; minimum slope is -80.03. 
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Evolution of max slope. n=4096, tfin=1.246, max|u_x|=80.03 



• Figure 7 Evolution of the maximum of |u x | for Run 1. 
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Final computed profile, n=4096, nu=0.1, max|u_x|=80.03, tfin=1.246 
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Log-linear evolution of spectral amplitudes, n=4096, nu=0.1 , t=0 to 1 .246 
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• Figure 9 Evolution of spectra] amplitudes. Log-linear plots corresponding to Run 
1. The filter level is set at 10 -12 . 
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Power law evolution, n=4096, nu=0.1 



Evolution of logarithmic decrement, n=4096, nu=0.1 



• Figure 11 (a) Evolution of the power law estimate of the spectrum; (b) evolution 
of the logarithmic decrement. Both for Run 1. 
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Evolution of the "enstrophy", n=4096, nu=0.1 , tfin=1 .246 



• Figure 12 Evolution of the enstrophy, u*(x , t)dx , corresponding to Run 1. 
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u(x;t) 



• Figure 13 Run 2. Evolution of u(x,t): Negative diffusion case v = —0.1; initial 
condition obtained by integrating uo(x) = sin(x) backwards in time; profiles are 
plotted every 0.05 time units and the final profile at t = 0.91 is also included. 
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Log-linear evolution of spectral amplitudes, n=4096, nu=-0.1 , t=0->0.91 
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